### Analysis for NAFTA hypotheticals
### Generates Table 1 in manuscript
### Last Revised and commented: 9 June 2016, JBS
##########################################################################################

# clear workspace
rm(list=ls())

# set working directory (needs to be changed by user); path will depend on where the folder is stored within in the file structure
setwd("~/Dropbox/IO imperfect monitoring/GLS_II_REP/")

# set random number seed for exact replication
set.seed(1234)

# load libraries
library(msm)

# load key functions for running the agent-based model
source("functions.R")

# generate a space with 100 units
x <- rep(0,100)

# set the number of simulations to 300
sims <- 300

# create holding bins for results from the simulations of the four scenarios
res1 <- vector("list",sims)
res2 <- vector("list",sims)
res3 <- vector("list",sims)
res4 <- vector("list",sims)
res5 <- vector("list",sims)
res6 <- vector("list",sims)


# Define positions and configurations; Mex=73, Can=89, US=90, Chile=79
# Rescaled as discussed in paper.
config1 <- c(46,78,80) # US/Canada/Mexico
config2 <- c(46,58,78,80) #US Canada Mexico Chile

# for each scenario, run the enlarge function
for (a in 1:sims){	
	### US/CAN/MEX high misperception
	res1[[a]] <-enlarge(space = x, members= config1, sdnmem = 6, a = 0.5, d=0.5, sdp = 2, R = 150)	
	### US/CAN/MEX low misperception
	res2[[a]] <-enlarge(space = x, members= config1, sdnmem = 6, a = 0.5, d=0.5, sdp = 4, R = 150)	
	### US/CAN/MEX/CHILE high misperception
	res3[[a]] <-enlarge(space = x, members= config2, sdnmem = 6, a = 0.67, d=0.5, sdp = 2, R = 150)
	### US/CAN/MEX/CHILE low misperception
	res4[[a]] <-enlarge(space = x, members= config2, sdnmem = 6, a = 0.67, d=0.5, sdp = 4, R = 150)
	### US/CAN/MEX high misperception, higher admit threshold
	res5[[a]] <-enlarge(space = x, members= config1, sdnmem = 6, a = 0.67, d=0.5, sdp = 2, R = 150)
	### US/CAN/MEX high misperception, higher admit threshold
	res6[[a]] <-enlarge(space = x, members= config1, sdnmem = 6, a = 0.67, d=0.5, sdp = 4, R = 150)
}

# Define function to calculate change variables
change <- function(sims,res,config){
	
	out <- matrix(NA,sims,3)
	colnames(out) <-c("newmems","medshift","newmempos")

	for (i in 1:sims){
		out[i,1]<-length(res[[i]])-length(config)
		out[i,2] <- median(res[[i]])-median(config)
		out[i,3] <- mean(res[[i]][-match(config,res[[i]])])
	}	
	out
}

# Apply change function to results
changeres1 <- change(sims,res1,config1) 
changeres2 <- change(sims,res2,config1) 
changeres3 <- change(sims,res3,config2) 
changeres4 <- change(sims,res4,config2) 
changeres5 <- change(sims,res5,config1)
changeres6 <- change(sims,res6,config1)

# Calculate means for changes
# These are reported in Table 1 of the paper
colMeans(changeres1,na.rm=T) # Table 1 upper right
colMeans(changeres2,na.rm=T) # Table 1 upper left
colMeans(changeres3,na.rm=T) # Table 1 lower right
colMeans(changeres4,na.rm=T) # Table 1 lower left

## Results for higher admission threshold in US/MEX/CAN agreement
colMeans(changeres5,na.rm=T)
colMeans(changeres6,na.rm=T)
